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Abstract 

During cell migration, cells become polarized, change their shape, 
and move in response to various internal and external cues. Cell 
polarization is defined through the spatio-temporal organization of 
molecules such as PI3K or small GTPases, and is determined by in¬ 
tracellular signaling networks. It results in directional forces through 
actin polymerization and myosin contractions. Many existing mathe¬ 
matical models of cell polarization are formulated in terms of reaction- 
diffusion systems of interacting molecules, and are often defined in one 
or two spatial dimensions. In this paper, we introduce a 3D reaction- 
diffusion model of interacting molecules in a single cell, and find that 
cell geometry has an important role affecting the capability of a cell to 
polarize, or change polarization when an external signal changes direc¬ 
tion. Our results suggest a geometrical argument why more roundish 
cells can repolarize more effectively than cells which are elongated 
along the direction of the original stimulus, and thus enable roundish 
cells to turn faster, as has been observed in experiments. On the 
other hand, elongated cells preferentially polarize along their main 
axis even when a gradient stimulus appears from another direction. 
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Furthermore, our 3D model can accurately capture the effect of bind¬ 
ing and unbinding of important regulators of cell polarization to and 
from the cell membrane. This spatial separation of membrane and 
cytosol, not possible to capture in ID or 2D models, leads to marked 
differences of our model from comparable lower-dimensional models. 


1 Introduction 

The ability to migrate is one of the fundamental properties of cells and is 
observed in both single-celled organisms as well as multicellular organisms in 
development, tissue maintenance, and in disease progression. For effective, 
directional migration, cells need to have the capability to sense and respond 
to various migratory signals, such as bacteria reacting to nutrients or other 
attractants or repellents BEI. cells being guided to their correct location 
during embryonic development |3] or immune cells migrating towards loca¬ 
tions of injury or infection [1]. Furthermore, cell migration plays a prominent 
role in diseases such as cancer, where the majority of deaths are caused by 
metastases. Hence migration, invasion and metastasis are considered dehning 
properties of cancer nu- 

External stimuli affecting cell migration include biochemical signals |7j 
or mechanical interactions with the environment [5113 HU]. One particularly 
interesting feature is the capability of many cells to detect spatial variations 
in the concentrations of biochemicals and to migrate towards, or away from 
the sources of such chemicals. Often, the gradients of those chemoattrac¬ 
tants or repellents have a small slope, so cells need a mechanism to detect 
and magnify external biochemical stimuli m- Such gradient detection then 
enables the cells to develop a polarized state with a well dehned front and 
back. To this purpose, chemical signals need to be translated into the gen¬ 
eration of mechanical forces [12] , which ultimately enable the cell to migrate 
in the direction dehned by the polarized state. 

In the last few decades, researchers have discovered and studied a large 
number of key molecules understood to play an important role in the sensing 
of chemical stimuli as well as the subsequent polarization, regulation of the 
actin cytoskeleton and generation of mechanical forces [13] • Among these 
molecules are small GTPases [m ITS] . PI3K, PTEN, phosphatidylinositols 
(PIPs), [13 HZl US] , Arp2/3 [13120] and Cohlin [2T1122] . 

To understand the complexity of those pathways of interacting molecules. 
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as well as to understand the mechanisms of sensing external gradients and 
polarizing a cell, a large number of mathematical models of gradient sensing 
and cell polarization have been developed (see [23l|23] for reviews). Whereas 
some of these mathematical models try to explain the general principles of 
signal detection, amplihcation and polarization [23 [221 EZ] , others attempt 
to explicitly model the dynamics and interactions between some of the most 
important involved molecules [2E1 EHl EHl EH EH EH EH E5] • Many of these 
models are formulated in terms of reaction-diffusion partial differential equa¬ 
tions (PDEs) and make use of ideas such as pattern formation, which have 
been applied to biology for many years [2Sl EH EHl IM]- Alternative modeling 
approaches to cell polarization include |1D], where thermodynamic consid¬ 
erations were used to predict polarization, im, where the effects of the in¬ 
terplay of biochemistry and mechanics on polarization were investigated, or 
|42j . where stochastic cell polarization was considered. The majority of these 
mathematical models have been formulated, or at least tested, in one or two 
spatial dimensions. Simulating a model in those lower dimensions greatly de¬ 
creases the computational costs, and might seem justihed if one is modeling 
cell migration on 2D substrates or in quasi-one-dimensional scenarios such 
as the detection of a ID chemical gradient. 

However, when the cell has an irregular shape, it is not a priori clear 
that a lower dimensional mathematical model can be used. Furthermore, in 
3D in vitro experiments or in vivo, stimuli can appear from all directions. 
An additional complication is the spatial organization of the key molecules 
behind cell polarization and migration: some of the regulators of the actin 
cytoskeleton, like the Arp2/3 complex, are soluble in the cytosol, whereas 
others such as phosphatidylinositols, are bound to the membrane. Moreover, 
some molecules such as the small GTPases can be both membrane bound and 
soluble, and this binding is influenced by the presence of other regulators such 
as guanine dissociation inhibitors (GDIs) [13]- Some mathematical models 
such as [H] have studied the influence of cell geometries on cell polarization 
and migration in two dimensions, the role of cell shape on signaling |1H , the 
mechanical effects of shape on cell migration [IH], the effect of cell shape on 
stress hber polarization m , or the effect of signaling on cell shape [Ml 129] , see 
also the review [50] . A mathematical model focusing specihcally on the effect 
of 3D shape on cell polarization, taking into account a whole polarization 
pathway, has, to our knowledge, not been investigated. 

In this paper, we are studying the effect of the cell shape on gradient 
sensing and cell polarization in a 3D mathematical reaction-diffusion model 
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of key molecules involved in polarization. In section we introduce a 3D 
model of GTPase molecules binding and unbinding from the membrane. We 
then use this framework to generalize an earlier pathway model investigating 
the dynamics of the small GTPases Rac, Rho and Gdc42, as well as PIP, PIP2 
and PIPS, in a ID context |33], and show how our 3D model can be reduced 
to such a ID model. Then, in section we explore how geometry affects 
the polarization capability and timescales of polarization. We first compare 
results of the 3D model with simple rectangular geometries to the ID limit 
and highlight similarities and key differences. We also find that varying 
the membrane binding and unbinding rate can change the cell polarization 
behavior. These rates are altered by the presence of GDI molecules, and 
while the dynamics of GDI molecules is not included in the present model, 
our results suggest how GDI molecules will affect cell polarization. We then 
show how cells with the same volume and length can have vastly different 
polarization behavior if they have different geometries. Finally, since in vivo 
migratory stimuli rarely appear constant in time and space, but dynamically 
change directions and strength, we investigate how cells react to changes 
in stimulus, and how this reaction is influenced by geometry. We find that 
if ellipsoidal cells are initially polarized along their main axis, they cannot 
adapt to a new stimulus perpendicular to their main axis as efficiently as 
symmetric, roundish cells. This gives a purely geometrical explanation of the 
fact that roundish, amoeboid cells can quickly turn and adopt to new stimuli. 
Furthermore, the ellipsoid cells preferentially polarize along their main axis 
even if the stimulus gradient is not aligned with this axis. The results in 
this paper thus predict that cell shape is an important factor influencing the 
ability of a cell to sense external signals, polarize and ultimately migrate. 


2 Models 

In this section we are introducing a 3D model of cell polarization, and discuss 
the relation to analogous ID models. In section 2.1, we define a model which 
consistently describes the binding and unbinding of a molecule to and from 
the membrane. Then, section uses this membrane-cytosol interaction 
model for inactive GTPases and includes activation of the membrane-bound 
GTPases, interactions of the three important small GTPases Rho, Rac and 
Gdc42 as well as interaction with phosphatidylinositols. Finally, in section 
2.31 we discuss how to reduce our 3D model to a ID model. 
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2.1 3D Membrane-Cytosol Interaction Model 

We denote by Gc the density of molecules which are freely diffusing in the 
cytosol, measured in moles per volume, and by Gm the density of membrane- 
bound molecules, measured in moles per area. The unbinding rate from the 
membrane is denoted by koff-, and the effective binding rate is konLj, where 
kon is a conventional rate with dimensions of inverse time, and Lj is the length 
scale dehning the region of the cytosol adjacent to the membrane which is 
accessible to the membrane-binding reaction. The diffusion coefficients for 
diffusion in the cytosol or on the membrane, respectively, are denoted Dc 
and Dm- Then, G^ and Gc evolve according to the following PDEs: 


dt 

aGc(fc,t) 

dt 

DqCji^ vGcij'my t'j 


D + konLiGJyTmG) koffGni(^Tj-fi, 


DcVlG,{n,t) 

konk'lGcij'myt^ kof fGmij'mG') ■ ( 1 ) 


The boundary condition for Gc ensures conservation of the number of molecules 
under binding and unbinding, and Bn is the unit outward normal vector at the 
membrane, so is the projection of the gradient on the normal vector. 

fm and fc denote points on the membrane or in the cytosol, respectively, and 
Vy, V| denote the volume and surface Laplace operators (otherwise known 
as Laplace-Beltrami operator, or Laplacian), respectively. Similar models 
as ([^ have been used in ^1] in the context of diffusion-driven instabilities. 
The boundary condition is also similar to the boundary conditions chosen 
to model the flux through a membrane, as done, for instance, in [52]. Such 
boundary conditions are known as Kedem-Katchalsky boundary conditions. 
A more detailed discussion and derivation of those equations is provided in 
the section SI of the supplementary information. We note that from ([^, 
it follows that when G^ and Gc are in equilibrium, and are homogeneously 
distributed, then the fraction / of membrane-bound molecules is given by 


/ = 


kn 


u \ U y ’ 

^on \ ^of f r.c 


( 2 ) 


''Off Li S 

where V is the volume and S the surface area of the cell. 
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Figure 1: The molecular pathway considered in this paper incorporates the 
GTPases Rac, Rho and Cdc42, and the phosphatidylinositols PIP, PIP2 and 
PIPS. The subscripts, a, mi and ci denote active, membrane-bound inactive 
and cytosolic inactive GTPases, whereas the PIPs are all membrane-bound. 
We consider the interactions of active GTPases and PIPs, the activation and 
deactivation of membrane-bound GTPases, and the binding and unbinding 
from the membrane of the inactive GTPases. 
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2.2 Pathway Model 

We now focus on a 3D cell polarization model incorporating the three GT- 
Pases Rac, Rho and Cdc42, as well as the three phosphatidylinositols PIPl, 
PIP2 and PIPS, as dynamic quantities. We use the molecular interactions as 
shown in Fig. which were previously considered in a ID model |33]- PIP, 
PIP2 and PIPS are all assumed to be purely membrane bound, whereas the 
three GTPases Rho, Rac and Gdc42 are assumed to exist in active and inac¬ 
tive membrane-bound forms, indicated by subscripts a and mi, respectively, 
as well as an inactive soluble form, which can diffuse in the cytosol and which 
is indicated by a subscript ci. The membrane binding and unbinding of the 
inactive forms is described as outlined in section 12.11 The full model is thus 
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described by equations (§. 
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Here, the activation functions Iq for the three GTPases are given by 
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( 4 ) 


and 6g are the deactivation rates. The signal SRad^m, t) is dehned to sim¬ 
ulate the effect of membrane receptor stimulation of Rac, i.e. it increases 
the Rac activation rate in a spatial way, and is thus dehned on points fm on 
the membrane. Typically, we will choose a function monotonously increas¬ 
ing along the direction of an external growth factor stimulus. More details 
about the relation of this model to the ID model of [33] are found in the 
supplementary information, section S3 In Table STl we also include a full 


list of the parameters appearing in our model dehned by equations ® and 

0 - 


2.3 ID Reduction 


We now consider the reduction of equation Q to a cylindrical cell of length 
L and radius R, where we assume cylindrical symmetry and no strong spatial 
dependence in the radial direction of the cylinder. Then, equation Q reduces 
to 

= DudlG^izR) + 2kon^G,{z,t) - koffGm{^,t), 

^ DcdjGciz, t) + (^koffGm{z, t) - 2kon^Gc{z, t)^ . (5) 


Here, Gm and Gc are the densities in one spatial dimension obtained from re¬ 
ducing Gm and Gc by Gm{.z, t) ~ 27rRGm(0, z, t), Gc{z, t) ^ 7ri?^Gc(r, 0, z, t), 
using cylindrical coordinates with radius r, angle 0 and axis ; 2 . In the deriva¬ 
tion, we have made use of the assumptions 


RLrkn 


< D, 


c, 


koffR^ 

6 


< D, 


C-, 


( 6 ) 
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which indicate that radial diffusion is fast, so concentrations equilibrate fast 
in the radial direction. Note that while equations ([^ are defined on a ID 
spatial domain defined by the length of the cylinder, so 0 < 2 ; < L, the 
cylinder radius is implicitly present in the sense that the membrane-binding 
rate kon is effectively renormalized by the inverse of the cylinder radius R. 
If we consider a cell with a given volume V = then, while maintain¬ 

ing a cylindrical shape, increasing the length L of the cylinder will result 
in a decrease of the radius R. Hence, the effective membrane binding rate 
‘^kon^ will increase. This makes intuitive sense as for a longer and thinner 
cylinder, proportionally more molecules in the cytosol are close to the mem¬ 
brane. Indeed, the thin layer of width Lj around the membrane, which is 
the region of the cytosol accessible to membrane-binding of the molecules, 
becomes larger for smaller R under fixed cylinder volume. The fraction / of 
membrane-bound molecules obtained for the ID cell is given by 


/ 


kn 


kon T k, 


off 


V 


47rLiL 


(7) 


which differs from the result of [33] • More details of the derivation of the 
ID limit are presented in the supplementary information, section S2.1[ and 
section S2.2 gives an analogous derivation for the reduction to two spatial 
dimensions. 


3 Results 


In section 3.1, we consider quasi-one-dimensional cells in our 3D framework 


and compare this with established ID models. Then, in section |3.2[ we will 
investigate the role of the membrane unbinding rate on cell polarization. 


Finally, in section 3.3 we will investigate how 3D geometry can influence 


the capability and timescales of cells to polarize and to repolarize when the 
external signal is changing directions. 


3.1 Polarization of a Quasi-One-Dimensional Cell 

To compare to the ID model [33], we are now investigating a scenario of 
emerging polarization where we start with initially homogeneous concentra¬ 
tions of all molecules, which are then perturbed by a large spike in active 
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(a) t = Os (b) t = 20s (c) t = 40s (d) t = 200s 


Figure 2: Transient polarization. We show the concentration of active Cdc42 
on the membrane for a rectangular cell with side lengths L,w,d = 20, 8,5/im. 
Active Rac is perturbed at the top of the cell, leading to a brief polarized 
state which then fades away with progressing time. 


Rac at one end of the long cell in a symmetric way depending only on the 
direction of the longest extent of the cell. First of all, we checked our code 
on a cuboid-shaped cell with side lengths L,w,d = 20, 8, 5fim, since a cuboid 
presents the most straight-forward generalization of a ID geometry. Fig. 
shows a typical time evolution of a GTPase concentration, here Cdc42, after 
the initial Rac stimulus is applied at time f = 0 at the top of the cell, which 
is then removed. We observe that at time t = 20s it looks as if the cell 
could polarize, but the strength of polarization fades away and is completely 
absent at time t = 200s. As we are interested in studying the effect of ge¬ 
ometry in this paper, we vary the length of the cell, fixing the cell volume to 
V = 800/im^. 

Fig-i shows this cell with different lengths, L = 40, 80/im. In each case 
we perturbed active Rac at the top of the cell and show the active Cdc42 
concentrations after t = 200s. We see that in both cases a stable polarization 
pattern is established. This is in contrast to the case of the shorter cell with 
L = 20/im, which, as shown in Fig. has no signs of polarization after 
t = 200s. Fig. I^in the supplementary information shows results from the 


same simulations but focuses on the time series of active Cdc42 at the front 
and back of the cell. These results are compatible with the observation in [33] 
that length can change the bifurcation behavior and increase the polarization 
sensitivity. However, our results are different for several reasons: we take the 
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(a) L = 40^m, t = 200s (b) L = 80^m, t = 200s 


Figure 3: Persistent polarization. The same cell and setup as in Fig. but 
elongated to L = 40 and 80/rm. We see that, contrary to the cell with L = 
20/rm shown in Fig. 2(d), after 200s a stable polarized state is maintained. 
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Figure 4: Length dependence of active Rac in homogeneous conditions (a) 


and in the presence of a gradient (b), in which case the concentrations at 
both the front and back are shown. In each case we compare the ID model 
from [33j with our 3D model, and the parameters are aligned so that, with 
homogeneous conditions and L = 20/im, both models coincide. 


hniteness of the membrane binding and unbinding rates into account; the 
fraction of membrane-bound GTPases which we derived in the section 12.31 
is different from the one used in |33]; |33] combined the inactive membrane- 
bound and cytosolic forms into one inactive form whereas we do not perform 
this approximation in the 3D model; [33] measured the membrane-bound 
particles in moles per volume, whereas we use moles per unit area, which is 
important as we fix the volume, but by changing the length also change the 
surface area of the cells. 

In Fig. |4(a)[ we show active Rac levels in a homogeneous, steady state 
setting, and in Fig. 4(b)| active Rac is shown in the presence of a constant 
linear gradient stimulus 200s after this stimulus is initially applied, for dif¬ 
ferent lengths of the cell in the ID and 3D models. As before, the volume 
of the cell is hxed at V = SOO/im^. We have chosen the parameters of the 
3D model such that at the base length of L = 20/im, we get agreement with 
the ID model and homogeneous conditions. Fig. 4(a) shows that the steady 


state values obtained in the homogeneous case differ signihcantly when the 
length is changed. Moreover, when the gradient is applied, the results pre¬ 
sented in Fig. 4(b) confirm that the differences in Rac active concentrations 
at the front and back of the cell can differ markedly between the ID and 3D 
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Figure 5: Geometry influences Polarization strength. As a measure of po¬ 
larization strength, we show the difference of active Rac on the membrane at 
the front and at the back of the cell, for cells of different rectangular shapes. 
Each plot shows four graphs for cells of different lengths, L = lO/im (orange), 
20/rm (gray), 30/rm (blue), 40/rm (green). In (a) the volume is hxed to be 
V = 800/rm^, and the sides w > d perpendicular to L of the rectangle are 

hxed in dependence on the asymmetry factor As as w = d = 

In (b), the volume dependence is shown while w = d = y 


model. In the scenario shown, in the 3D model, the difference of active Rac 
at the front and the back decreases with increasing length, whereas in the 
ID model it increases with increasing length. 

Furthermore, in a ID model one cannot accurately take into account 
variations in the directions perpendicular to the main axis. In Fig 


5 (a) 


we 

investigate the impact of the asymmetry factor As = y, which describes the 
asymmetry of the directions perpendicular to the length L, on polarization. 
We hnd that, for cells of different length, higher asymmetry decreases polar¬ 
ization strength, measured in terms of the difference of active Rac between 
the front and the back of the cell. Then, in Fig. 5(b), we investigate the 


volume dependence, and hnd that generally, increasing the volume V of the 
cell increases polarization strength for cells of different lengths. A main effect 
of changing either volume or asymmetry is that this will change the volume 
to surface ratio, which then affects effective activation and inactivation rates 
as well as effective membrane binding and unbinding rates. Furthermore, the 
effective diffusion rates are changed when volume or asymmetry change. 
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10 20 
kof,[l/s] 

(a) L = 20fim 

Figure 6: Membrane unbinding influences polarization strength. As a mea¬ 
sure of polarization strength, we compare the difference in concentrations 
of active Rac between the front and the back of the cell as a function of 
koff for different values of the Rac activation rate Iri and two lengths, 


L = 20/im and (b) L = 40/im, in all cases 200 seconds after an initial stimu¬ 
lus of active Rac is applied at the front of the cell. In (a), only the cell with 


Iri = 0.0005 shown in red is polarized, and all three other values lead to 
completely unpolarized states. 


3.2 Role of Membrane Unbinding Rates 

We now investigate the dependence of our model on an important new pa¬ 
rameter typically not considered in previous models, that is, the membrane 
unbinding rate koff. Its associated binding rate kon is fixed via relation (|^. 
It is of physiological importance, as GDI molecules mediate the sequestration 
of GTPases into the cytosol [^, and hence their dysregulation will change 
binding and unbinding rates. Hence, it is important to know how the model 
predictions change when these rates are varied. Fig. shows the difference of 
active Rac between the front and back, which is a measure of the polarization 
strength, of a rectangular cell as a function of koff 200s after an initial stim¬ 
ulus of Rac at the front. We show this dependence for different cell lengths 
and values of the Rac activation rate Iri, as these parameters where shown in 


to be important parameters affecting cell polarization. Fig. 6(a) shows 
results from a cell of length L = 20/rm, whereas Fig. 6(b) shows results 
from a cell of L = 40/im. In each case, we show plots for four different Rac 
activation rates Iri- First, we confirm qualitatively the observation of 
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that intermediate ranges of Iri can lead to a polarized states, or are more 
strongly polarized. Furthermore, the shorter cell does not polarize as easily 
as the longer cell, as in Fig. |6(a)| only the cell with Iri = 0.0005 is polarized. 
We also see that koff is positively associated with polarization strength, and 
the cells with very small values of koff = 0.1s“^ do not, or only weakly, po¬ 
larize. Note that the ID limiting case requires koff ^ ^ for 

L = 20,40/xm, respectively, so in neither case are the approximations applied 
in [33j necessarily expected to be accurate. This constraint is discussed along 


equation (25) in the supplementary information. We also note that for most. 


but not all parameters checked, the polarization strength saturates at koff 
values of the order of magnitude of ls“^. The results in Fig. [^conhrm that 
the membrane unbinding rate is an important parameter which can influence 
the capability of a cell to polarize. 

3.3 Influence of Geometry on Polarization 

We now investigate how cell shape influences the ability of the cell to polarize, 
lose polarization or repolarize when the direction of a signal changes in time. 
Here we include some effects which cannot be investigated with a ID model. 


3.3.1 Influence of Geometry on Initial Polarization 

In many experiments, cells present in vastly different shapes. For instance, 
as shown in Fig. a cancer cell is extravasating through a vascular lumen 
As it does, it transitions from a nearly spherical shape (Fig. 7(a)), 


into one consisting of a spherical region inside the lumen, spreading into 
a broad set of protrusions in the extracellular matrix outside of the lumen 
(Fig. |7(b) ). The two parts of the cell are connected by a thin neck-like region 
reaching through the endothelium, barely visible in this single confocal slice, 
but typically being of about 1 or 2/im in diameter. While a full model of 
the extravsation process would need to take complexities such as the change 
in environment from the blood-hlled inside of the lumen to the extracellular 
matrix outside of the lumen into account, our current model provides, with 
limitations, an understanding of what effect complex cell shapes similar to 
those seen in Figure [7b] would have on the polarization behavior of cells. 

We compare a cell with two different shapes: First as a single ellipsoid. 


Figs. 8(a) (c), then, as two ellipsoids joined by a thin neck 1.3/im in diam¬ 


eter between the ellipsoids. Figs. 8(d) (f)] For a better comparison we keep 
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(a) (b) 


Figure 7: Cells can appear in vastly different shapes. Here, a cancer cell 
(green) extravastating from inside a vessel of endothelial cells (purple) into 
surrounding extracelullar matrix (black) is shown, as observed in [5l] (Repro¬ 
duced by permission of The Royal Society of Chemistry). In (a) the cancer 
cell appears nearly spherical, while it is still fully inside the vessel lumen 
and has not started to extravasate. When it is in the process of extravasa¬ 
tion through the endothelium, it narrows dramatically at the endothelium, 
connected only through a thin neck region (b). Part of the cell remains in 
the lumen, but much of it has already spread outside of the lumen into the 
extracellular matrix. 


the length and volume of the two conhgurations the same, so that the main 
difference between the two cases is the thinning, and the spreading of one 
half of the cell, similar as seen in the extravasating cell outside of the lumen 
in Figure We see that for both shapes, the cell is polarizing at t = 5s in 
response to the stimulus. However, at t = 100s, the single-ellipsoid cell has 
lost its polarization. Fig. 8(c), whereas the extravasating cell maintains a 


strongly polarized state. Fig. 8 (f), such that active Cdc42 is mainly concen¬ 
trated in the part of the cell outside of the lumen. This could explain the 
formations of hlopodia, known to be directed by Cdc42, almost exclusively 
outside of the lumen. However, as mentioned before, the current model does 
not take all complexities during the extravasation process into account so 
further work is required to investigate if shape alone, or a combination with 
other effects such as the presence of ECM molecules outside of the lumen, 
are responsible for the observed behavior. 
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Figure 8: Active Cdc42 for cells of different shapes: an ellipsoid cell (top 
row, (a)-(c)), and a cell composed of two thinly connected ellipsoids (bottom 
row, (d)-(f)) at times of 5, 10, and 100s. A stimulus gradient is applied to the 
initially homogeneous cells. Both cells initially polarize at t = 5s, but the cell 
with the shape of a single ellipsoid loses its polarization at f = 100s, whereas 
the cell composed of two ellipsoids is able to maintain its polarization. 


3.3.2 Response of Cell to a Change in Stimulus Direction 

We now investigate how cell shape can influence the response of a cell to a 
change in the direction of a stimulus. This will allow us to make predictions 
with this model beyond the response to unidirectional stimuli, which exist 
typically in in vitro setups such as classical Boyden chambers [55] , or modern 
microfluidic platforms [5511571158] 159] . However, many cases of cell migration 
in vivo are more complex due to tissue heterogeneity as well as temporal and 
spatial changes in stimuli [60|. The scenario presented in Fig. |^goes towards 
an understanding of how cells respond to changes in stimuli, and how this 
response is affected by cell shape. Active Rac is shown for a cell of a volume 
oi V = 800pm^ for three different ellipsoidal conhgurations: With a main 
axis of 11.5/im (spherical, (a)-(d)), 15/im ((e)-(h)) and 20/im ((i)-(l)). Ini¬ 
tially a stable polarized state is obtained in all three conhgurations, which is 
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Figure 9: Active Rac on the membrane is shown at different times for the 
same cell with different shapes, where the Rac activation rate in the first 
100 s increases linearly along the long axis of the ellipsoid (from lower right 
corner to upper left corner), and from then on, it is rotated by 90 degrees 
and now increases linearly along a short axis of the ellipsoids (from the lower 
left corner to the upper right corner). In all cases, the volume of the ellipsoid 
cells is fixed as V = 800/im^, the main axis is 11.5/im (spherical, (a)-(d)), 
15/im ((e)-(h)) and 20/im ((i)-(l)), and the other two axes are of the same 
length. Comparing the different shapes, we see that only the spherical cell can 
completely polarize into the new stimulus direction, whereas the cells with 
ellipsoidal shapes will form a stable pattern which points into a direction in 
between the original and final stimulus direction. 
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Figure 10: As in Fig. active Rac on the membrane is shown for different 
times and cells of different shapes, but here, the Rac activation rate in the 
hrst 100s increases linearly along a short axis of the ellipsoid (from lower 
left corner to upper right corner), and from then on, it is rotated by 45 
degrees and increases linearly in between the long and a short axis of the 
ellipsoids (from bottom to top). In all cases, the volume of the ellipsoid cells 
is fixed as R = 800pm^, the main axis is 11.5/im (spherical, (a)-(d)), 15/im 
((e)-(h)) and 20pm ((i)-(l)), and the other two axes are of the same length. 
Comparing the different shapes, we see that the spherical cell can completely 
polarize into the new stimulus direction, whereas the cells with the elongated 
ellipsoidal shapes will rotate their internal polarization direction further than 
the stimulus direction, so that the polarization is more aligned with their long 
axis. 
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shown at t = 100s just before the activation gradient direction is changed. 
The spherical cell always maintains a main direction, but the direction of 
polarization rotates towards the new stimulus direction. At t = 180s the cell 
is mainly polarized into the new stimulus direction (c), and at t = 300s it 
is completely repolarized. Interestingly, the ellipsoidal conhgurations never 
reach a polarized state aligned with the new stimulus. Instead, the new 
polarization direction is somewhere in between the original and hnal stimu¬ 
lus direction, and the larger the ratio of the ellipsoid axes is, the closer the 
new polarization direction will remain to the original polarization direction. 
This would suggest that a more symmetric, spherical cell is able to repolar¬ 
ize faster and more efficiently than a long, thin cell, despite their chemical 
pathways being unaltered. In Fig. IT, we use the same setup as in Fig. 


but we initially polarize along one of the short ellipsoid axes and then rotate 
by 45 degrees towards the long axis. The spherical cell can, as before, adapt 
to the new stimulus direction. Now, the ellipsoidal cells can rotate faster, 
and, interestingly, they rotate their internal polarization state by more than 
45 degrees towards the long axis. Indeed, the longest ellipsoidal cell is, after 
t = 300s, almost fully polarized along the long axis (Fig. 10(1)). Supple¬ 
mentary Fig. S3 shows the same setup, where, like in Fig. f^we initially 


polarize along a short axis, but then fully rotate by 90 degrees towards the 
long axis. Here, all cells can repolarize towards the new direction, and the 
ellipsoidal cells repolarize faster than the spheroidal cells. It is well known in 
the literature that roundish, amoeboid cells can quickly polarize and adapt 
to new stimuli, in contrast to more elongated cells such as mesenchymal cells 
|61j . Futhermore, in 3D matrices, elongated mesenchymal cells can migrate 
persistently along the hbrous structures. Whereas traditionally mesenchymal 
and amoeboid cells denoted different cell types, such as hbroblasts on the one 
hand, and dictyostelium or neutrophils on the other hand, recent work has 
also focused on a switch between those migratory modes for the same cell 
types. Most of this work has focused on alterations of biochemical pathways 
to describe the switch between mesenchymal and amoeboid migrations. In 
|62j . the role of MMPs in this switch was investigated, in |H3] an important 
role of Rho and ROCK was discovered, and in [M] LIM kinase was implicated 
in this switch. In [65], cell aspect ratios were taken as the factor determin¬ 
ing whether a cell migrates in a mesenchymal or amoeboid way, with the 
mesenchymal cell being more persistent, and interstitial flow was shown to 
affects the switch between these migratory modes. Here, our argument shows 
that, without any changes in the biochemical pathways, roundish cells are 
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expected to repolarize towards a new stimulus more quickly than elongated 
cells purely because of their different shape. On the other hand, elongated 
cells preferentially migrate in the direction of their longest extent, even if the 
stimulus appears in a slightly different direction, giving them an increased 
persistence. 


4 Conclusion 

Motivated by the fact that cells, both in in vitro and in vivo environments, 
present with greatly varying shapes, in this paper we have investigated how 
cell shape influences gradient detection and cell polarization. For molecules 
such as small GTPases, which can exist both in a membrane bound form 
as well as sequestered in the cytosol, the ratio of cell volume to surface 
area can influence important properties such as the fraction of membrane- 
bound molecules, as seen in equation ([^. Similar arguments were provided 
in |in]. However, all GTPases are affected in this way, and they are partially 
inhibitory to each other, so it is not a priori clear how a full model with inter¬ 
acting molecules would behave under change of cell shape. Furthermore, we 
found that even in the dimensionally reduced models, the higher dimensions 
implicitly appeared in the reduced equation by renormalizing the membrane 
binding rate (equation ([^). We also found conditions on the parameters 
which, when satisfied, justify the use of the lower dimensional models (equa¬ 
tion (j^). These conditions are quite general and should be useful both to 
check the validity of older models as well as for the development of new mod¬ 
els which involve the interactions of molecules between the membrane and 
the cytosol. 

In our model, we have only considered constant binding and unbind¬ 
ing rates of a molecule to and from the membrane, and investigated how 
these rates affect the polarization behavior of a cell (Fig. |^. In real cells, 
the binding and unbinding of small GTPases to and from the membrane 
as well as their activation and deactivation are influenced by the presence 
of a large number of different molecules such as GDIs, guanine nucleotide 
exchange factors (GEFs) and GTPase activating proteins (GAPs) |53]. It 
would be interesting to include the effect of these regulating molecules in 
our model. However, at present, there is a lack of good quantitative data 
regarding the spatio-temporal regulation of these molecules, so we postpone 
such investigations for future research. What Fig. confirmed is that the 
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binding/unbinding rate can influence the ability of a cell to polarize, and as 
GDI molecules are expected to modify those rates the model predicts that 
the presence or absence of these molecules will also affect the polarization 
behavior. 

We then extended an established cell polarization pathway [3S], which 
was previously investigated in a ID model of HeLa cells, to our 3D model. 
The purpose of choosing [33j for comparison was that in this model, the effect 
of changing the cell length was implicitly taken into account via a modifi¬ 
cation of the fraction of membrane-bound GTPases, whereas most other ID 
models did not consider any geometric effects at all. Hence, a first test of 
our model was to reproduce some results of [33] and highlight quantitative 
and qualitative differences. Furthermore, as mentioned above, we are, to our 
knowledge, the first to explicitly consider the binding/unbinding dynamics 
of GTPases to the membrane and show the influence of those parameters on 
the polarization behavior of the cell in Fig. 

We then explored scenarios which the existing lower dimensional models 
could not capture. First, we compared the polarization behavior of two cells 
with the same volume and length, one ellipsoid, and one cell composed of 
two connected ellipsoids. Fig. The second shape was motivated by shapes 
observed during cancer cell extravasation [53] , where the cellular environment 
can impose different extreme shapes on the cells. This result is important 
whenever one is trying to compare theoretical results obtained from simplified 
lower dimensional models to experiments, indicating that one has to take cell 
shape into account. We expect that the behavior of pathways other than 
those describing polarization would also be affected by cell shape in a similar 
manner, if the principle mechanism of polarization is mediated by similar 
reaction-diffusion models as the ones used here. 

As long as the polarization stimulus is coming from only one direction, 
and provided the parameter constraints (|^ are satisfied, ID models could 
still be derived which take into account if the cross section of the cell along 
the stimulus direction is relatively constant. However, purely one-directional 
stimuli are idealistic and in vivo different stimuli can appear from all direc¬ 
tions and change in time. In the study of such effects we have to use 3D 
models such as ours. In Fig. [^we simulated a stimulus which was changing 
directions over time. We found that cell shape is an important factor which 
predicts how fast and strong a cell can adapt to the new direction of stim¬ 
ulus. Indeed, the spherical cell was able to change the internal polarization 
direction smoothly towards the new stimulus direction, whereas cells which 
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are elongated along the previous direction of stimulus preferentially stayed 
polarized in a direction close to the original stimulus. This is compatible 
with experimental hndings that roundish amoeboid cells are typically faster 
to adapt to new stimuli than mesenchymal cells, which are typically more 
elongated. An explicit test of the model prediction could be conducted, for 
instance, using a microfluidics platform where one can change the direction 
of an external growth factor gradient over time, and measure the response of 
some tagged internal molecule associated with polarization for varying cell 
shapes. This would be a step toward an understanding of cell polarization 
under temporally and spatially varying conditions as typically present in vivo 
|60j . From a theoretical point of view, it would be interesting to include the 
effect of dynamical changes of shape through coupling of mechanics with our 
biochemical pathways, as these dynamical changes have also been shown to 
affect polarization behavior on longer time-scales in a 2D model |44j . 

In summary, the results in this paper predict the importance of cell shape 
on polarization of cells, indicate in which cases the use of lower dimensional 
models is justihed, and demonstrate when a full 3D model such as ours needs 
to be used to model and predict experimental results. 
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Supplementary Information 


SI Derivation and Consistency of the Mem¬ 
brane Binding-Unbinding Model 


We model a molecule which is diffusing in a cell, can bind to and unbind 
from the cell membrane, and diffuse on the membrane when bound. Let 
Gc be the concentration of this molecule in the cytosol, i.e. the inside of a 
cell. We denote the domain of the cytosol by V C which is a smooth 
Riemannian manifold with the metric induced from the Euclidean metric in 
Likewise, Gm is the concentration of the same molecule when bound 
to the membrane, which is defined as the boundary oi V, S = dV, and is 
an orientable Riemannian manifold. The membrane-bound molecules can 
unbind, and the molecules in the cytosol can bind to the membrane, with 
rates kon and koff, and where Lj is a length scale associated with the binding 
range of a sequestered molecule to the membrane. Furthermore, Dm and Dq 
denote the diffusion coefficients for diffusion on the membrane and in the 
cytosol, respectively. Our equations are given by ([^, which we repeat here 
for convenience: 


dGmij'my t) 
dt 

dG^{fc, t) 
dt 

Gm(?’m) 0) 
Gc(rc, 0) 


^cV^Ge(fe,t), 

konL iG (JyTmit^ koffGjYi(j^iYi,t^, 

Glirn.), 

G°(re). (8) 


Here, V|, Vy denote the Laplace operators (otherwise denoted as Laplacian, 
or Laplace-Beltrami operator) on S and V, respectively, and are defined in 
the usual way on Riemannian manifolds [21! • Furthermore, e„ denotes the 
uniquely defined unit outwards normal vector on the surface, and fc G R, 
fm G S. Hence, e„VyGc(rm,t) denotes the projection of the gradient of Gc 
on the unit normal vector on the surface. We have imposed the outwards 
normal flux in such a way that it matches the binding and unbinding reactions 
and preserves total particle numbers. Furthermore, G^{fm),G^{fc) denote 
functions defining the initial conditions, and naturally G^ does not need any 
boundary conditions, as it is defined on a snrface without boundary. 
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Sl.l Particle Number Conservation 

We now show that the equations given in (|^ conserve the number of particles. 
The total amount of molecules is given by 


N= GcdV+ / GmdS. 
Jv Js 


( 9 ) 


This total amount of molecules is conserved by choice of boundary condition: 


dN 


dG, 


dV + 


dGr 


-dS 


J dt J dt 
= / Dod^v,(gra4yGMV + j(D,Mvs{gradsG^) + k^L.G^ - K„G^)dS 

= j Dcen{gradvGc)dS + j{DMdivs{gradsGm) + KnLiGc-koffGm)dS 

= j DcCnigradvGcjdS + J{0 + konLjGc — koffG 

m )dS 

/ ( kQrfiL jGq T T kQrfiL jGq kfy^ ^G^g^dS 


= 0 


( 10 ) 


In the hrst line of the derivation, we have simply plugged in the time deriva¬ 
tives of Gc and Gm from equation ([^. Here, we have used that (on Rie- 
mannian manifold, independent of the coordinate system) we can write the 
Laplace operator as divergence of a gradient, where the subindices indicate 
the corresponding manifold in which divergence or gradient are calculated. 
We do not need the precise dehnition of gradient, divergence or Laplacian on 
those manifold, we only need the fact that the divergence theorem applies. 
Indeed, we apply the divergence theorem in the second step, changing from 
an integral of a divergence of the gradient of Gc over the whole cell to an inte¬ 
gral of the normal flux over the boundary. Then, in the third step, we apply 
the divergence theorem to the divergence of the gradient of Gm- However, 
the surface does not have a boundary, so the divergence theorem immediately 
gives zero for this term. In the hnal step, we plug in the boundary condition 
from equation ([^ for the normal flux of Gc at the boundary, and obtain our 
hnal result, that the total amount of particles is conserved. 
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SI.2 Global Invariance of Boundary Condition 

Global conservation of particles is still ensured by adding a Laplacian of Gm 
to the flux boundary conditions: 

( 11 ) 

It is immediately clear that this preserves total particle numbers for any A, 
for if we integrate the boundary condition over the whole boundary, this 
term drops out by application of the divergence theorem over a manifold 
with empty boundary. However, for local conservation, we should not keep 
this term: If we choose kon = ^o// = 0, then m does not bind or unbind from 
the membrane at all and is conserved on its own (and not just the sum of 
bound and unbound molecules). However, keeping boundary condition ( [II| 
with A 7 ^ 0 would result in a flux of c even in that case. On physical grounds, 
we have to impose Neumann no-flux boundary conditions and the right-hand 
side of ( pH| should be zero. Hence, only A = 0 ensures local conservation of 
particles. 


SI.3 Derivation from a Model with Finite Binding Ra¬ 
dius 

We can consider a generalization of equations 0.il where the change of the 
membrane density m{fm,t) is affected by all molecules in the cytosol within 
a hnite radius Lj from the point fm- Then, the membrane binding-unbinding 
model is described by the equations 


Wt 


= DMVlG^{r^,t) 

+ ~k on / Gc{rc,t) - 

J\fm-rc\<Li 


dGc.{rcG) 

dt 


= DcVlG,{f,,t), 


■J \rm-rc\<Li 


We have assumed that all molecules which are with a distance of Lj to a 
point on the membrane fm are equally likely to be bound with a rate of 
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kon, which could be generalized further by including a kernel in the inte¬ 
gral such that molecules closer to the membrane are more likely to bind. 
However, if we assume that Lj is small (see the discussion on parame¬ 
ters in the supplementary information S3.2) so that Gc does not signih- 


cantly vary on this length scale, and that the membrane is not signihcantly 
curved on this scale, we can reduce the integral terms in (12). The in¬ 
tegral will thus be over a half-sphere with radius Lj, so we can simplify 
kon = kon^^'n'LjGc{f^,t) = konLiGc{fm,t), wliere we 

have identihed kon = ^T^Ljkon- Hence, (12) reduces to equations ([^, ([^. 


SI.4 Alternative Derivation from a Discrete Model 

We now give an alternative derivation of equation ([^ from a discrete model. 
We consider a small section of a cell near the cell membrane, so small that 
we can ignore the curvature of the membrane. Such section is shown in a 
schematic drawing on the right panel of Fig. SI, where the cell membrane 
is highlighted by the red surface. We are interested in the dynamics of the 
binding and unbinding of molecules to the membrane. Let Lj be the inter¬ 
action length such that when a molecule in the cytosol is within a distance 
less or equal to Lj of the membrane, there is a probability of binding this 
molecule to the membrane. The associated binding rate is denoted by kon- 
Likewise, unbinding is denoted koff. 

We dehne the domain of interest to be a cube of length L ^ Li^ which 
we discretize into equally spaced small cubes of size 5. Initially, we identify 
5 = Lj. Each cube is labeled by integer-valued indices {k,l,p), and the 
membrane is located at the boundary p = 0. Then M{k,l) denotes the 
number of membrane-bound molecules at the membrane segment adjacent 
to cube {k,l,0), and C{k,l,p) denotes the number of cytosolic molecules in 
the cube {k,l,p). We consider the following processes: In the inner part of 
the cytosol, unbound molecules can diffuse only. At the cube adjacent to 
the membrane, they can diffuse in parallel to the membrane or away from 
the membrane, or they can bind to the membrane. On the other hand, 
membrane-bound molecules can unbind, or diffuse on the membrane. The 
following equations describe the rate of changes of the average number of 







(a) Diffusion in cytosol 


(b) Binding to the membrane 



(c) Diffusion on membrane (d) Unbinding from the membrane 

Figure SI: Molecules in the cytosol can move from their current cube to any 
of the nearest-neighbor cubes (figure (a)), or, if they are within a distance of 
the interaction range Lj to the membrane, they can bind to the membrane 
(hgure (b)). Membrane-bound molecules can diffuse on the membrane only 
(£gure(c)), or unbind from the membrane (hgure (d)). 
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molecules: 


^ ^M{k + 1,0 + M{k - 1, 0 + M{k, / + 1) + M{k, / - 1) - 4M(fc, /)) + 

ot 0^ 

konC{k, I, 0) - koffM{k, 1), 

^ +1) Op) + C{k — 1, Op) + C{k,i + i,p) + C{k,i — i,p) 

+ C{k, l,p+l) + C{k, l,p-l) - 6C{k, /)), p > 0, 

= ^{C{k + 1,1,0) + C{k - 1,1,0) + C{k, I+ 1,0) + C{k, I-1,0) 
ot 6^ 

+ C{k, 1, 1) - 5C{k, /)) + koffM{k, 0 - konC{k, I, 0). (13) 


We now want to study the continuum limit of those equations. First, we add 
an artificial layer of cubes at p = — 1, such that 

Cik, I, -1) := ^ {koffM{k, 1) - konC{k, I, 0)) + C{k, I, 0) . (14) 

The benefit of this layer is that we can now combine the last two equations 
of (13) into 

^ +1)Op) + c{k — 1,Op) + c{k,i + i,p) -\- c{k,i — i,p) 

+ C{k,l,p+l) + C{k,l,p-l) -0C{k,l)), p>-l, (15) 


which is supplemented by (14). From those equations, it is straight-forward 
to take the continuum limit 5 —)■ 0. We define Gc(a;, y, z) = to be the 

density of unbound molecules at a point [x, y, z), where we identify {x, y, z) = 
{6k, 61, 6p), and likewise Gm{x, y) = is the surface density of molecules 

bound to the membrane. Then, the continuum limit 5 —?■ 0 gives the following 
equation: 

^G'"(^,^)-= T)^V|Gm(a:,i/) + fcon [ Gc{x,y,z)dz - koffGm{x,y) 


D. 


dt 

dGcjx, y,z) 
dt 

dG c{x,y,0) 
dz 


c- 


= DcV^yG,{x,y) 

I'Ll 

= Kn / Gc(x, y, z)dz - koffGmix, y ). 


(16) 


30 














We have introduced the 2 and 3D surface or volume Laplace operators, V| 
and Vy, respectively. If we assume that the concentration in the cytosol 
does not vary much on a length scale of L/, we can reduce these equations 
to the local model given by equation ([^. 


SI.5 Binding-Unbinding Equilibrium 


We now consider the equilibrium condition between the binding and un¬ 
binding to and from the membrane, which is obtained from equation ([^ by 
equating LikonGc = koffGm- The rates of binding and unbinding are set by 
koff and ^ respectively, where R is the cylinder radius. If / denotes the 
fraction of membrane-bound to total concentration. 


/ = 


JsGmdS 


f^GJV + f^GmdS 

then, in a homogeneous equilibrium, we get 

k 

j, rvnn 


h -L h ^ ’ 

^on 1“ 1^0ff ijs 


as 


argued in the main text in equation ([^. 


(17) 


( 18 ) 


S2 Dimensional Reduction of the 3D Model 

We are now deriving the dimensional reduction of the 3D membrane bind¬ 
ing/unbinding model Q. 

S2.1 ID Reduction 

We consider a cylindrical cell where the height of the cell is L and the radius 
is R, such that L ^ R. We choose cylindrical coordinates (r, 0, z) such that 
r G [0,i?], z G [0, L] and cylindrical symmetry, so our helds and Gc do 
not depend on 0 G [0,27r]. As Gm is the concentration of membrane bound 
molecules, it is only dehned at the boundary of the cell located at r = R and 
z = 0, L. Furthermore, the dependence of Gc on r is weak, relative to the 
dependence on z, due to fast radial diffusion due to L S> i?. We completely 
neglect the r dependence of G^, as this would only matter at a; = 0, L, and 
would be weak, similar to the weakness of the r dependence of Gc. 
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The PDEs in cylindrical coordinates then become 


dC (z f) 

^ = DMdlGUz,t) + konLlG,{R,Z,t) - koffGm{z,t) 

^ Gc(r,2:,t) 

Rc^rGzkQfiLjGc(^R^ z T kQjjGfyii^zR'j^ 

DcdzGc{r, 0 , t) = konLiGc{r, 0 , t) - koffGm{0, t), 

DcdzGc{r, L, t) = —konLiGdr, L, t) + koffGm{L-, t)- (19) 


Then, let us dehne the following ID densities 


r-i? /■27r 


Gc{z,t) = 


Gc(r, z, t)rd(j)dr 


'0 ^0 
rR 


= 2ti 


Gc(r, t)rdr, 


i>2tt 


Gm{z,t) = / Gm{z,t)Rd(l) = 27lRGm{z,t). 


( 20 ) 


As R is assumed to be small such that diffusion in the radial direction is faster 
than other timescales in the problem, radial diffusion will quickly homogenize 
Gc in the radial direction even in the case when the initial conditions have a 
strong r dependence. Hence, ignoring potential fast transient changes of Gc, 
we consider Gc to have a weak dependence on r and expand Gc{z, t) to get 

Gc(2:, t) = 2 ti f {Gc{R, z, t) + drGc{R, z,t){r — R) + ...) rdr 

Jo 

/ d2 d3 1 

2tt ( Gc{R,zR)— - — — {koffGm{z,t) - konLiGdR, z,t)) 

V / D JJc 


= 7rRUl + 


RLi 

3Dc 


kon ) Gc(.R, Z, f) 


koffR^ 

6Dc 


Gm{z^ f). 


( 21 ) 


Here, we have expanded about r = R and used the boundary condition. If we 
further assume that the radius R and the interaction range Lj are small such 
that the equation ([^ holds, then we can approximate further (see section 


S3.2 for a discussion of those parameters) 

Gc{z^ t) ^ TrR‘^Gc{R, z, t). 


( 22 ) 
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The ID densities follow the differential eqnations 


^ = DMdlGm{z,t) + 2nRkonLiG,{R,z,t)) - k^ffG^zR), 

= DMd^Gm{z, t) + 2kon^Gc{z, t) - koffGm{z, t), 

JrC 

(^dlGRz, ^) + (^drrdrGRr, z, t)^ 2Trrdr^ , 

= Dc (dlGc{z, t) + 2T:RdrG,{R, z, f)) , 

konLlGc^yZ^t^^ 

GRzR^ , 


= Dc (^lGc{zR) + 2 tiR^ {koffm{z,t) - 
= DcdlGc{z,t) + (^koffGm{z,t) - 2koJ^ 


1 


(23) 


which were given in the main text in eqnation (|^. Mass conservation is 
ensnred by accompanying those eqnations by Nenmann no-flnx bonndary 
conditions. Note that while these eqnations are perfectly ID PDEs, with 
the spatial domain dehned by the length of the cylinder, the cylinder radins 
is still felt in the sense that the membrane-binding rate kon is effectively 
renormalized by the inverse of the cylinder radins R. 


S2.1.1 Equilibrium in One Spatial Dimensional 


We consider the binding-nnbinding eqnilibrinm condition for the cylindrical 
cell in the ID limit, as done in section S1.5| for a general 3D cell. We obtain 
for the fraction / of membrane-bonnd molecnles 


/ = 


kn 


T ffoT^ 


h//2Lr kon + k, 


off 


V 

i-KLjL 


(24) 


as argned in eqnation Q. Notice that the dependence of / on D and L is 
different to the one fonnd in [33], and we have checked that this difference 
is not becanse of the use of rectangular, rather than cylindrical cells. The 
binding/unbinding equilibrium and simultaneous ID limit are valid if the 


33 











parameters satisfy the constraints 




1 Dc 

> koff > 

‘^KnLj ^ Dc 
R ^ T^' 


(25) 


These constraints ensnre diffnsion in the radial direction occnrs on the fastest 
timescale to ensnre there is only a weak radial dependence, and a ID limit 
is justihed. On the other hand, binding and nnbinding are faster than diffn¬ 
sion along L so the binding/nnbinding eqnilibrinm is jnstihed. We can then 
introdnce an effective diffnsion coefficient Dmc = IDm + (1 ~ f)Dc, similar 
as in |29ll33], so that the total concentration = Gc{z,t) -|- Grn{z,t) 

simply evolves by the standard ID diffnsion eqnation 


z,t) ^ d^G^°\zR) 

dz^ 

We have thns rednced the original system of two conpled PDEs in three 
spatial dimensions to a single PDE in one spatial dimension. Note that 
other potential reactions in the system wonld need to be modihed by / ac¬ 
cordingly, and, if they are present, the timescales of associated with those 
reactions need to be compared to the timescales of radial diffnsion and bind¬ 
ing/nnbinding to jnstify the rednction of the complete model with reactions 
to lower dimensions. 




S2.2 2D Reduction 


Let ns consider a flat cell, which, for simplicity, we take to be a disk of radins 
R and height h. Hence, it is natnral to choose cylindrical coordinates, and 
the fnll model is described by (19), with L replaced by h. We now consider 
the limit h R. Then, we can rewrite the 2 ;-dependence of the Laplacian as 

d‘^Gc{r,(j),z,t) 5^Gc(r,0, h) - 9^Gc(r,0,O) 

R>c -~ R>r 
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kpff 

h 

konL I 

h 


h 

(Gm(r, h,f) -F Gm(r,0,f)) 
(Gc(r,0, h,t) + Gc(r, 0,0, f)). 


(27) 
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If the concentrations only depend very weakly on z, then we can simplify the 
system by introducing 


Gc(r, 0, f) = / Gc(r, 0, z, t)dzdr = hGc(r, 0, t), 


Gm(r, 0,f) = 2Gm(r, 0,t), 


(28) 


and get 


= Z^ArV^G^(r,0,t) + 2fc„„^a(r,0,t) - koffG^{r,cl),t), 

(9G (r d) f) L 

—= DcVpGc{r, 0, t) - 2konj^Gc{r, 0, t) + koffGjn{r, 0, t). (29) 


Here, Vp denotes the conventional 2D Laplace operator in polar coordinates. 
Similarly to the ID case, the scaling of the parameters is different. Similarly 
to the ID case, we hnd that the parameters of the reduced geometry, in 
this case, the cylinder height h, renormalize the effective membrane-binding 
coefficient. As before, we consider the steady-state solution where / denotes 
the fraction of membrane-bound molecules. Hence, for the oblate cylinder 
we get 


/ 


+ ^off2L, 


(30) 


S3 3D Polarization Pathway 


In this section, we discuss in detail how our 3D pathway model discussed in 
is obtained and relates to the ID model discussed in 


section 


2.2 


S3.1 Model Setup 

The ID model of [ 33 ] was motivated by experiments where cells were con¬ 
strained in effective ID geometries. It was assumed that the approximate 
3D geometry is rectangular with length scales L ^ w > d, with an initial 
length of L = 20/im. The volume of a cell in the experimental paper M, 
which uses the model of [ 33 ], was given as approximately V = 800/um^. The 
supplementary information of (33] mentions d = 0.2/um, which seems a bit 
small for a real cell and would also imply that w = 200/um at the given 
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volume. We hence compare to a cell with base measure L,w,d = 20, 8, S/im. 
In [67], the cells were about 80/im long (see e.g. hgure S4 in [67|), which 
would, at the same volume, be compatible with L,w,d = 80, 5, 2/im. Such 


cell seems also compatible with the ID limit as described in the section |2.3 
of the main text. In [331EZ] the change in cell length was taken into account 
by changing the fraction of membrane bound to unbound inactive molecules 
via an equation similar, but slightly different, to our equation ([^. The 3D 
model will automatically take into account any geometry change. The ba¬ 
sic equations for the evolution of the three GTPases now takes into account 
that the each GTPase can exist in an active, membrane bound form Gma, a 
membrane bound, inactive form G^i and a form Gc which diffuses inactively 
through the cytosol. These follow the principal scheme 
FtC 

+ — Og'^ 


dt 

dGmi 

dt 

—DcenVyGc 


= DmV|g, 


= DM^sGrm - IcGr 


= DcV^Ge, 

konLlGc kgj j^G-uii- 


T f^cGma T konLjGc kgffG 


of 


(31) 


Here, Iq represents the activation, and 5g the deactivation rate, whereas 
komkoff denote the binding and unbinding rates as in the section]^ Note 
that equations (31) are slightly different from the equations given in the 
appendix of [33], which were used to motivate the ID model from a 3D 
perspective. 

We note that to account for the proper localization of the membrane 
bound and unbound species G^a, Gmi and Gc, we measure Gc in Molar, 
but Gma, Gmi In Whereas this latter measure is not often chosen in 

experiments, as usually total cell concentrations are measured, this is never¬ 
theless the physically more meaningful measure, as Gma, Gmi denote number 
molecules per two dimensional membrane area, and this choice ensures that 
our equations and the dimensional reductions have the correct units. To 
compare with [33], we will hence multiply the concentrations of active GT¬ 
Pases with ^ = 4/3/im, the fraction of volume to surface area for 

the above mentioned rectangular cell of basic length L,w,d = 20,8,5/im. 
With this setup, it is straight-forward to generalize the ID model from [33] . 
summarized in hgure in our 3D context, and one obtains equations (|^. 
All coefficients apart from the membrane binding and unbinding rates are 
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taken from IB7] , but those which multiply a membrane density are mul¬ 
tiplied by the factor Furthermore, to compare the activation rates given 
with ours we need to divide them with the fraction of bound to 


m 


total inactive molecules /, as we separately consider bound and unbound 
inactive GTPases. We now summarize the use of the spatially dependent 
Rac stimulus Sjiac which appears in the activation rate for Rac, equation 
(|^ in most simulations shown in section If the x-axis denotes the direc¬ 
tion of the stimulus, then typically we assumed a linear stimulus of strength 


SRac — 0.5Iri 


in Table SI 


20^m' 


where Iri is the baseline Rac activation rate as given 


The baseline length of 20/im is chosen to match results from 
earlier works, as described above. In hgures and ^ such gradient was 
used but then rotated towards the indicated axis at f = 100s. 


S3.2 Parameter Estimation 

We now investigate when the use of the ID model is justihed for the case 
where the molecule is a small GTPase. We have Dq = lOOpm^ /s [291 E7j. 
Furthermore, we estimate that the interaction range of the binding reaction, 
L/, approximately corresponds to the size of the molecules. We have a molec¬ 
ular weight of the small GTPases of about 21kDa. Exact size determination 
of proteins is tricky [HH], but here we only need a rough estimate, which 
gives that we have a volume of R = = 3 * = 30nm^. 

Hence, the interaction length scale is on the order of a few nanometers. We 
put Lj = 2nm, and this estimate is similar to stimations made in similar 
contexts [69]. Furthermore, we can safely assume that R 3> T/ for real¬ 
istic cell geometries. Then, unless kon S> koff, of the two requirements 
<C 1, 1, the hrst one automatically holds provided the sec¬ 

ond one does. From I7D1 we can estimate that koff should be faster than 
koff = 0.06s“^, as the combined deactivation/unbinding rate (this combined 
rate is denoted koff in [70|) is of this magnitude. However, the actual bind¬ 
ing and unbinding rates are influenced by the presence of other regulators 
such as GDI molecules [13] and might be different for GTP and GDP bound 
GTPases, and is hence also influenced by the presence of GEFs and GAPs. 
Here, we focus on rough estimates and use the above numbers to derive a 
limit for the radial length scale of 
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As long as koff is not too large this condition is satisfied for realistic cell 
dimensions. However, if ko/f should be significantly larger than estimated 
above then this limit might be hard to satisfy. 

Now we look at the steady-state assumption between bound and unbound 
GTPase. For this, we have 


With koff = 0.06s“^, we get L ^ 30/rm would safely satisfy this constraint. 
However, most likely koff is significantly larger so the steady-state assump¬ 
tion is most likely valid for shorter cells as well. 

All other parameters used in equations (|^ are summarized in Table SI 


Lacking accurate measurements of koff-, it is commonly believed that koff is 
much faster larger than the deactivation rate [TOl |33]. As the deactivation 
rates were estimated in [33] to be we take koff = 10s“h While the 

total unbinding/deactivation rate was estimated in [7Q] to be much smaller 
than ls“^, we stick to those values here as we first would like to compare our 
3D model to the ID model of |33|. In the subsection 3.2 in the main text 


we study the influence of varying koff on the polarization behavior of the 
cell. We also need to determine the combination of parameters konLi- From 
the estimates of Dc-, Dm, combined with the estimate that at baseline length 
of L = 20fim the diffusion coefficient for total inactive GTPases (bound 
and unbound) is 50fim‘^/s, we get that about half of the inactive GTPases 
molecules are typically membrane bound. This equilibrium value can then 
be used to deduce konLi via equation ([^. 


S3.3 Implementation of the 3D Model 

All simulations of the 3D model were performed in GOMSOL Multiphysics 
5.1 (GOMSOL, Inc, Burlington, MA) using the General PDF model frame¬ 
work. We remark that the default solver occasionally produced too large 
time steps, requiring us to manually limit the maximal time step depending 
on the model parameters. It is also necessary to choose a fine mesh for good 
spatial resolution in several cases, for instance, when the cytosolic species 
vary sharply at the membrane. In most cases, the predefined Mesh Element 
settings ’Finer’ or ’Extra Fine’ were sufficient to ensure spatial convergence. 
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Dimensional parameter 

Dc 

Dm 

Dp 

RhOtot,2 

Cdctot,2 


Estimate 
100 /im^ s“^ 
1/rm^ s“^ 
5/rm^ s“^ 
lOnmol m“^ 
4nmol m“^ 
3.4nmol m“^ 


Psb 

0.2nmol m ^ 

Sr 

Is-i 

Sp 

Is-' 

Sc 

Is-' 

Spi 

0.21s-' 

^21 

0.021s-' 

h 

'^max 

2.8s-' 

kp2 

2.1s-' 

a 

1.3s-' 

Pp 

0.011s-' 

G 

0.03s-' 

Iri 

0.4pM s-' 

Ir2 

0.4pM s-' 

Ir\io 

13.2/iMs-' 

Icdc 

5.9pMs-' 

IPi 

14nmol m-^ s 

(l\ 

1.7nmol m-^ 

Oj2 

1.3nmol m-^ 

fcpiSK 

0.084s-' 

fcpi3K 

0.00072s-' 

/CpTEN 

0.432s-' 

fl 

1 

^off 

10s-' 

konLi 

13.3pm s-' 


Table SI: Parameters of the cell polarization model equation (|^, inferred 
for a cell with dimensions L,w,d = 20, 8, bfim, so all values which multiply 
membrane concentrations are rescaled by the factor ^ = |/im. This means 
Furthermore, a, J/ji, 1r 2 ^Rho and Icdc are multiplied by 
the baseline fraction of bound inactive molecules /. We have also rounded 
the parameters as appropriate. 
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(a) L = 20fim 



(b) L = ‘iOfim 



Figure S2: The time evolution of concentration of active Cdc42 in time at 
the front and back of the rectangular cells as shown in hgures and 
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(a) t=100s 

(e) t=100s 

(i) t=100s 



(f) t=140s 

(j) t=140s 



(c) t=180s 

% 

(g) t=180s 

% 

(k) t=180s 



(1) t=300s 



Figure S3: As in figures active Rac on the membrane is shown for 

different times and cells of different shapes, but here, the Rac activation 
rate in the hrst 100s increases linearly along a short axis of the ellipsoid 
(from lower left corner to upper right corner), and from then on, it increases 
linearly along the long axis of the ellipsoids (from the lower right corner to 
the upper left corner). In all cases, the volume of the ellipsoid cells is hxed 
as V = 800/rm^, the main axis is 11.5/im (spherical, (a)-(d)), 15/im ((e)-(h)) 
and 20/rm ((i)-(l)), and the other two axes are of the same length. Cells of all 
shapes are able to adapt to their new stimulus direction, with the elongated 
cells being slightly faster. 
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